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Abstract 

This survey paper assesses the status of compressible Euler and Navier-Stokes solvers on unstruc- 
tured grids. Different spatial and temporal discretization options for steady and unsteady flows are 
discussed. The integration of these components into an overall framework to solve practical problems 
is addressed. Issues such as grid adaptation, higher order methods, hybrid discretizations and parallel 
computing are briefly discussed. Finally, some outstanding issues and future research directions are 
presented. 
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1 Introduction 

Computational Fluid Dynamics (CFD) has evolved rapidly as a discipline and is increasingly being used 
to complement the wind tunnel, especially in preliminary design [50]. Based largely on the mathemat- 
ical foundations laid among others by Lax [88] and Godunov [52], the field has come into its own in 
the last decade. Great advances have been made in the areas of spatial discretization, grid generation 
and solution strategies. Tremendous advances in computer architecture and networking speeds have 
contributed to the field significantly as well. By the mid 80's, many different groups around the world 
were able to compute three-dimensional flows over simple realistic aerodynamic configurations. The 
grids employed were of the body-fitted or structured grid type. One-dimensional models were extended 
to deal with multiple dimensions in a natural way because of the structure by using the so-called gener- 
alized coordinates [161]. However, the task of generating structured grids about complex configurations 
presented a serious challenge. The widely-used multi-block structured grid approach solves this problem 
by tessellating the domain between the body and the far-field into simple logically rectangular blocks, 
so that structured grids can be generated easily within each block [173, 81]. The automation of the 
blocking and the grid generation process are difficult tasks that are continually being refined. Another 
powerful approach uses overlapping or chimera grids [19, 110]. Here, structured grids generated about 
the different components, are allowed to overlap. Automation of blocking, grid generation and the 
preprocessing required for deriving the interpolation operators are continually being improved. 

The desire to compute flows over complex configurations also spawned a surge of activity in the 
area of unstructured grids. The term “unstructured grids” will be used primarily in this paper to mean 
grids composed of simplices , which are triangles in two dimensions and tetrahedra in three dimensions. 
Unstructured grids have ahvays been used in finite element circles, but have become popular in the 
finite volume community only fairly recently. They provide flexibility for tessellating about complex 
geometries and for adapting to flow features, such as shocks and boundary layers. An underlying 
premise is that unstructured grid generation is far more automatable than are the tasks associated with 
multi-block structured grid generation. 

It should be mentioned that the flow' solver only constitutes a part of the overall solution method- 
ology. Other tasks, such as grid generation, interfacing with geometry packages, visualization and 
post-processing are equally important, and demand considerable skills and resources. The outline of 
this survey paper on unstructured grid flow' solvers is as follow's. Section 2 presents the governing 
equations in integral form. Section 3 reviews some popular finite volume spatial discretization schemes. 
Section 4 review's the finite element approach for spatial discretization. Section 5 review's the turbulence 
models in vogue for unstructured grid computations. Section 6 presents the various alternatives avail- 
able for time-discretization for steady and unsteady problems. Section 7 addresses the important topic 
of grid adaptation. Sections 8 and 9 examine higher order accurate schemes and hybrid discretizations. 
Section 10 briefly examines parallel computing issues. The paper concludes wdth a listing of issues to be 
resolved and some future research directions. In addition to the papers cited in this article, the reader is 
encouraged to read, especially regarding spatial discretization, the AGARD Report [1] on unstructured 
grid methods for advection dominated flow's. 

2 Governing equations 

The equations governing compressible fluid flow T in integral form for a control volume V(/) with boundary 
S(t) are given by 

4- [ Wdv + l [ F(W, to, s) — G(W, VW, n)]da = 0, ( 1 ) 
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In the formulas given above p is the density, V is the velocity vector with Cartesian components V,, e is 
the specific total energy, n is the outward unit normal vector of the boundary S(t) and s is the velocity 
vector of the boundary. Also, ft is the molecular viscosity, A is the bulk viscosity related to p by Stokes’ 
hypothesis, A = — 2/3/i, 7 is the identity tensor, T is the stress tensor and D is the deformation tensor 
given by 
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where VI j denotes the partial derivative of the ith component of V with respect to the Cartesian 
coordinate Xj, i.e. = ||C 0 stands for the divergence of V given by V l}l with the usual summation 
convention, q is the heat flux given by Fourier’s law 


q = — A'VT, (3) 

where K is the thermal conductivity of the fluid and T is the temperature. These equations are 
augmented by the equation of state, which for a perfect gas is given by 

P=(7-l ){pe~\p\V\ 2 ) (4) 

Eqn. (1) represents the conservation laws for the mass, momentum (the Navier-Stokes equations) and 
energy. It holds for any volume and in particular, holds for a specific volume associated with each grid 
point, termed the control volume. 


3 Finite volume spatial discretization 

It is assumed that a grid about the geometry of interest has been generated by some suitable method. 
The various grid generation techniques in use are reviewed in the survey papers by Thompson and 
Weatherill [166] and Mavriplis [106]. On a given grid, one has at least two choices as to where to locate 
the variables, giving rise to the cell- vertex and the cell-centered approaches. In the cell-vertex approach, 
the variables are stored at the vertices of the grid, whereas in the cell-centered approach they are stored 
at the centroids of the cells. There is yet another approach that stores only the averages associated 
with control volumes. If the scheme possesses second order spatial accuracy or less, this approach is no 
different than a cell-centered approach; the higher order scheme is covered briefly in Section 8. 

The concept of using arbitrary control volumes to solve numerically the conservation laws was es- 
tablished by the late 70’s, at least in theory [89]. Jameson and Mavriplis [73] reported some of the 
earliest, results from solving the two-dimensional Euler equations on regular triangular grids that were 
obtained by subdividing quadrilateral grids. In a cell-centered setting, they extended much of what was 
established for finite volume schemes on structured grids [74, 67] to triangular grids. This included the 
constructions of central-difference-like approximation for the convective terms and a blend of dissipative 
terms to suppress odd-even decoupling and to capture discontinuities, and the incorporation of multi- 
grid ideas. Solutions and convergence rates of comparable quality to structured grids were obtained. 
Second order accuracy was demonstrated by using a sequence of uniform grids. In 1986, Jameson et 
al. [72] presented their paper dealing with inviscid transonic flow over a complete aircraft. The pa- 
per’s contributions included grid generation for complex geometries using the Delaunay triangulation 
approach and the development of a cell-vertex flow solver. The control volume for each vertex in the 
tetrahedral grid was taken to be the union of all tetrahedra sharing that vertex. They also showed 
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the correspondence between this and a Galerkin finite element procedure with piecewise-linear basis 
functions. The nondissipative approximation was augmented with dissipative terms. In addition, the 
paper also presented a first discussion on a scheme with positive coefficients for triangles, although the 
proposed scheme was only first order accurate. Figure 1 shows the surface Mach contours for inviscid 
transonic flow over a Boeing 747-200 taken from the paper of Jameson et al. [72]. The flow conditions 
are M 0 0 = 0.8 and a = 2.73°. The mesh used for the calculation contained 12,038 nodes (57914 tetra- 
hedra) and the calculation took about 25 minutes on a Cray XMP computer. Since this seminal effort, 
tremendous advances have been made in grid generation and vectorization of the flow solver and are 
outlined in the sequel [71]. 

In 1973, Boris and Book [25] introduced the Flux Corrected Transport (FCT) scheme that converted 
a first order accurate monotone scheme to a higher order scheme by adding limited amounts of anti- 
diffusive flux to prevent spurious oscillations. Harten [54] developed a mathematical formulation that 
gave rise to the Total Variation Diminishing (TVD) schemes which also added limited amounts of anti- 
diffusive fluxes to a first order scheme to achieve monotonicity. Van Leer [169] independently devised the 
concept of limiting when designing monotonicity preserving schemes for conservation laws. He devised 
the Monotonic Upstream-centered Scheme for Conservation Laws (MUSCL) schemes which relied on a 
piecewise-polynomial reconstruction procedure that enforced monotonicity principles by using nonlinear 
functions called limiters. The jumps at the interfaces were resolved by using a Riemann solver, much like 
Godunov’s scheme. Since the Riemann problem for the Euler equations requires an iterative procedure 
and is expensive, the quest was on for approximate Riemann solvers. Roe [142, 143] and Osher [120] 
constructed two such schemes that proved to be particularly good in concert with the MUSCL approach 
on structured grids. For an excellent discussion on the development of upwind schemes, see the text by 
Leveque [91]. By the mid 80k, MUSCL schemes were being used by a number of groups to compute 
aerodynamic flows on structured grids [164, 160]. These schemes were gradually being extended to deal 
with unstructured grids. 

Desideri and Dervieux [40] devised cell-vertex finite volume schemes for unstructured grids using 
MUSCL ideas. Given pointwise values at the cell vertices of a triangulation, they employed a recon- 
struction procedure that made use of gradients in neighboring triangles. The flux was constructed using 
Osher’s approximate Riemann solver. Limiters were used in a one-dimensional fashion but a multi- 
dimensional monotonicity principle was not satisfied. Nevertheless, respectable results were obtained 
for many inviscid computations. Lohner et al. [95] tested a FEM-FCT scheme for Euler and Navier- 
Stokes equations. Fezoui and Stoufflet [45] proposed and tested a class of implicit upwind schemes 
that utilized various upwind approximations. Whitaker et al. [185] constructed a similar scheme using 
Roe’s flux-difference splitting scheme and obtained results for many transonic and supersonic flow r s. 
Venkatakrishnan and Barth [176] constructed an upwind scheme for cell-centered triangular grids that 
also employed MUSCL ideas as did Batina [16], Knight [84] and Frink [47], However, none of these 
efforts really guaranteed the absence of oscillations in the multi-dimensional case. 

Barth and Jespersen [15] made a radical departure from one-dimensional thinking for satisfying 
monotonicity principles. They enunciated a monotonicity principle in multiple dimensions similar to 
that employed by Van Leer [169] and Spekreijse for structured grids [160], namely that the reconstructed 
distribution in the control volume be bounded everywhere by the values of the neighbors (including 
the vertex representing the control volume) and satisfied the principle by constructing a truly multi- 
dimensional limiter. They employed Roe’s approximate Riemann solver for the evolution phase. As 
an ultimate test, they computed oscillation-free solutions for transonic flow^ over an airfoil on a highly 
irregular mesh. They obtained solutions of comparable accuracy wfith both cell-centered and cell- vertex 
schemes. This MUSCL approach has been adopted since by other researchers [5, 7, 51]. The multi- 
dimensional limiter in [15] may be thought of as a generalization of the min-mod limiter and as such, 
leads to convergence difficulties. Venkatakrishnan [174] has analyzed this problem and has proposed 
modifications that ameliorate the situation at the expense of monotonicity. Aftosmis et al. [3] have 
found that the modifications proposed in [174] significantly improve the convergence as w^ell as solution 
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accuracy in test problems. Shu [153] has proposed the Total Variation Bounded (TVB) schemes which 
weaken the TVD property by allowing for small-scale oscillations although the motivation here is to 
avoid clipping smooth extrema. The effect of TVB modification on convergence has not been studied 
yet. 

Frink et al. [47, 49] developed an upwind cell-centered three-dimensional flow solver. They employed 
a weighted averaging procedure that interpolated variables from the centers of tetrahedra to the vertices 
and used these vertex values to compute the gradients within each tetrahedron. These gradients were 
used to interpolate the variables to the centers of the faces of the tetrahedra. This was followed by the 
use of Roe’s approximate Riemann solver. The reconstruction procedure was linear and therefore, not 
monotonicity preserving. Nevertheless, the averaging procedure seemed to introduce enough dissipation 
that good transonic flow solutions were obtained without using limiters. The reconstruction procedure 
has since been modified to be linearity-preserving [48] based on the work of Holmes and Connell [61]. 
Laminar viscous capability has been been added to this code. Figure 2a depicts the surface grid and 
“oil flow” pattern for laminar flow over a delta wing at - 0.3, a = 20.5°, Rei = 0.9 x 10 6 . 
The flow has been computed by Frink using the methodology described in [48]. The grid containing 
730,454 tetrahedral cells was generated by the Advancing- Layers Method [132] which produces thin- 
layer tetrahedral grids suitable for computing viscous flows. Evidence of primary, secondary and tertiary 
vortices may be seen in Figure 2a. Figure 2b shows the pressure profiles at four chord stations and the 
comparisons with the results from the structured grid code CFL3D [164] and the experimental data of 
Hummel [66]. 

Barth [11, 12] presented extensions of the schemes of [15] to three-dimensional inviscid and viscous 
flow computations. The main contributions of these papers included the discretization of the viscous 
terms using a finite element procedure and the use of edge-based data structures for discretizing the 
inviscid and viscous terms. In addition, least-squares and data-weighted procedures for the construction 
of gradients were outlined. Mavriplis [102, 105, 101] developed an explicit vertex-based finite element 
multigrid scheme for the two- and three-dimensional Euler and Navier-Stokes equations. He also pre- 
sented edge-based data structures. In [11, 105, 98, 97] a case is made for using cell- vertex approximations 
and edge-based data structures in three dimensions by examining the computational complexity and 
storage. 

The issue of cell-vertex vs cell-centered approximations is still an open one, particularly in three 
dimensions. In two dimensions the ratio of the number of cells to the number of vertices is 2 whereas 
in three dimensions this ratio could be arbitrarily large, although it is typically between 5 and 6 for 
nice tetrahedralizations. In the case of tetrahedral grids, the flux computations in a cell- vertex scheme 
can be cast as loops over edges w T hereas in a cell-centered scheme they are loops over triangular faces. 
The ratio of the number of faces to the number of edges in a typical tetrahedral grid is roughly 2. 
Therefore, it. is true, as argued in [11, 105] that on a given grid, the cell-centered approximation incurs 
considerably more computational effort and memory requirements. However, there is some evidence 
that on a given grid, the solution quality' of a cell-centered scheme is superior to that of a cell-vertex 
scheme [134]. This is more than likely due to the fact that for triangular/tetrahedral grids the control 
volumes in a cell-centered scheme are smaller than those in a cell-vertex scheme. It is not clear whether 
the cell- vertex scheme requires a grid that has as many vertices as the number of tetrahedra employed 
by a cell-centered scheme to achieve the same level of accuracy. If this is the case, the cell-centered 
scheme wall prevail because it requires much smaller grids to be generated. However, it appears that a 
cell-vertex scheme is better suited to computing the viscous fluxes, especially when the triangulations 
become highly irregular. As discussed later in this section, the viscous terms in a cell-vertex scheme 
are typically discretized using a Galerkin finite element approach, whereas a finite-volume viewpoint 
is adopted in a cell-centered scheme. It has been shown [11, 90] that a finite element discretization 
of a diffusion operator with linear basis functions obeys the discrete maximum principle if a Delaunay 
triangulation is used. This result is true only in two dimensions and does not hold in three dimensions. 
On the other hand, the finite volume discretization of the diffusion operator in a cell- centered setting 
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does not guarantee a discrete maximum principle in two or three dimensions. It is interesting to note 
that the centroidal control volumes that are used with cell-vertex schemes are not necessarily convex 
and can assume rather odd shapes, especially in the case of highly stretched grids. One alternative is 
to use the Voronoi cells, which are guaranteed to be convex [86], but this approach requires care at the 
boundaries. 

Significant developments have also taken place in the area of positive schemes. These require that 
the semi-discrete scheme be expressible in the form 


c ij( u j - 

ai jeVi 


( 5 ) 


where the Cij ' s are non-negative, and A r i denotes the set of neighbors of i. As observed in [72], if the 
Qj's are constants, the scheme is linear and hence only first order accurate. Higher order versions 
of such schemes have been developed [42, 33, 69, 70]. These may be regarded as generalizations of 
TVD schemes to conservation laws in multi-dimensions. Durlofsky et al. [42] have developed a cell- 
centered scheme that makes use of gradients in auxiliary triangles similar to the scheme of [176], but 
incorporates a multi dimensional limiter to prevent oscillations. Both the schemes have the drawback 
of requiring nice triangulations, since the auxiliary gradients could be ill-defined on highly nonuniform 
meshes. Jameson [69] has developed a theory for Local Extremum Diminishing (LED) and Essentially 
Local Extremum Diminishing (ELED) schemes leading to Symmetric and Upstream Limited Positive 
(SLIP and L T SLIP) schemes. He has also presented an analysis for schemes to resolve shocks with one 
interior point [70]. A simple scheme that meets the conditions is termed a Convective LTpwind and Split 
Pressure (CUSP) scheme, which is closely related to the AUSM scheme of Steffen and Liou [92]. The 
schemes have been tested on structured grids for inviscid and viscous flows [162]. 

A different avenue to realizing upwinding on unstructured grids is provided by the fluctuation- 
splitting schemes [39, 121]. In multi-dimensions, the weak link in the MUSCL formulation is the use of 
a directionally-split Riemann solver. This approach is routinely adopted, but can misinterpret features 
not aligned with the grid. To overcome this problem, various multi-dimensional Riemann solvers have 
been proposed in the literature [38, 125, 148, 146], but have not gained widespread acceptance. Rather 
than adopt the approach of reconstruction followed by an approximate Riemann solver to bring in 
the upwinding, fluctuation-splitting schemes consider the average time evolution of a complete cell (a 
triangle or a tetrahedron) with the unknowns located at its vertices, and then update the values at 
the vertices by the effect of linear wave solutions evolving the piecewise linear data over a cell. If the 
distribution formulas are unbiased, this leads to the well-known Ni's scheme [119], a centered scheme 
which requires the addition of dissipative terms. Several uwpind fluctuation-splitting schemes for scalar 
advection equation have been developed; these are compared in [121]. Extension of fluctuation-splitting 
schemes to systems of equations relies on the decomposition of flux divergence into scalar waves. The 
decomposition is not unique, and several strategies such as characteristic decomposition [38] and simple 
wave models [144, 125] have been investigated. Some of the advantages of the fluctuation-splitting 
approach are that they result in compact stencils and that they do not use ad-hoc dimensionally-split 
Riemann solvers. Some of the disadvantages include the fact that the scheme is only second order 
accurate at steady state, and that they bind the user to using a decomposition technique. Recently, a 
positive linearity-preserving fluctuation-splitting scheme has been constructed for the multi-dimensional 
scalar advection equation using limiters [155]. It differs from the other approaches in that a purely 
algebraic (not geometric) viewpoint is adopted. This approach has enabled Sidilkover [154] to devise a 
fluctuation-splitting scheme for the Euler equations that may be viewed as a genuine multi-dimensional 
extension of Roe’s scheme without having to resort to decomposition into scalar waves. A unique 
property of this scheme is that a nonlinear Gauss Seidel procedure can be used as a relaxation scheme. 

In finite volume schemes, a spatial (possibly monotonicity-preserving) approximation for the inviscid 
terms is combined with a centered approximation for the viscous fluxes. A basis for this “operator- 
splitting” approach is furnished by Roe[145] who showed by means of elegant analysis in one dimension 
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how a monotonicity-preserving discretization for the inviscid terms when combined with a centered 
approximation for the viscous terms suppresses spurious modes and also yields a nonoscillatory solution 
for all cell Reynolds numbers. On unstructured grids, the viscous fluxes are discretized by using either 
finite volume or finite element approaches. In the case of cell-centered discretizations Frink et al. [48] 
compute the first derivatives at the vertices of the triangulation, which are then averaged to obtain the 
viscous fluxes at the faces. In the case of cell-vertex schemes, the viscous terms are usually derived by 
adopting a Galerkin finite element approach [147, 11, 101], although there is an equivalent finite volume 
interpretation. This formulation results in a compact stencil involving only the nearest neighbors, and 
also exactly reproduces the gradient of a linear function. A simpler formulation has been proposed in 
[98]. Here, the viscous fluxes are computed by averaging the cell-vertex gradients, which are available 
from the reconstruction procedure for the inviscid terms. However, this method results in a larger 
stencil for the viscous terms compared to the finite element discretization and could be more diffusive. 


4 Finite element discretizations 

Finite element method (FEM) merits a separate discussion since its evolution in CFD has taken place in 
parallel with finite volume schemes. The rigorous treatment adopted in FEM has important implications 
for global error estimation and grid adaptation. Many researchers have shown the equivalence between 
finite volume and finite element approaches when piecewise-linear approximations are employed for the 
solution vector and the fluxes [147, 72, 11, 150]. 

The most powerful method in finite elements is the Galerkin method. Here, the solution is first 
expanded in a set of basis functions and the residual is made orthogonal to a set of test functions. W hen 
the basis and the test functions are the same, the method is termed a Galerkin method; otherwise it 
is called a Petrov-Galerkin method. The standard Galerkin method leads to a centered scheme and is 
unconditionally unstable for hyperbolic problems when combined with forward Euler discretization in 
time. Therefore, artificial viscosity lias to be added in some form to stabilize the procedure. 

In contrast with the finite volume methods, finite element practitioners have always preferred cell- 
vertex schemes since the global function is usually expressed as a summation of trial functions multiplied 
by the values at the vertices. These trial functions are typically assumed to be 1 at a vertex and zero at 
all other vertices. Bristeau et al. [28] computed many low Reynolds number compressible flows using a 
finite element method. Angrand and Dervieux [8] investigated several first and second order accurate 
explicit schemes. They added Lapidus-type artificial viscosity to a Galerkin method. Donea [41] derived 
the important Taylor-Galerkin family of schemes for the linear advection equation. Adopting an FEM 
approach, he showed how a Galerkin scheme (a centered scheme) could be stabilized by using a Taylor- 
series expansion for the time derivative . similar to the procedure used to derive the Lax-Wendroff 
scheme. He demonstrated that the resulting schemes had good phase error and dissipation properties 
and that they could be easily extended to multi-dimensions. Lohner et al. [96] developed and tested a 
two-step Taylor-Galerkin finite element method for the solution of the Euler equations. Since the use 
of artificial viscosity spread the discontinuities over several cells, they also employed an adaptive grid 
algorithm using local refinement. 

W r hile there are finite difference/volume equivalents to most of the finite element techniques des- 
cribed in the previous paragraph, a very important class of schemes has been derived in the finite element 
community that has no obvious counterpart in the finite difference/ volume approach. Included in this 
class are the Streamwise Upwind Petrov Galerkin (SUPG) or the Streamwise Diffusion (SD), and the 
Galerkin Least Squares methods. SUPG, devised originally by Hughes and Brooks [63] for the steady 
scalar advection- diffusion equation, and Galerkin least squares [64] methods have subsequently been ex- 
tended to deal with the compressible and incompressible Navier-Stokes equations [62], To the Galerkin 
discretization, they add terms proportional to the residual (unsteady residual for time-dependent prob- 
lems) to introduce dissipation. Unlike conventional artificial viscosity, this form of artificial viscosity 
does not compromise the order of accuracy of the scheme since the exact solution yields zero dissipation. 
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The global error in £2 norm can been shown to vary as h v+l / 2 for the advection-dominated case and as 
/i p+I for the diffusion-dominated case, where p is the degree of piecewise- polynomial basis function [62], 
Another method that has been extensively analyzed is the discontinuous Galerkin method [77, 33, 22]. 
Johnson et al. [78] have analyzed the SUPG scheme and have extended it to deal with time-dependent 
problems by adopting a discontinuous Galerkin procedure in time. A review and analysis of these and 
other finite element methods may be found in the text by Johnson [76]. All these methods are linear and 
are not suitable for capturing discontinuities. Hughes and Mallet [65] have proposed a shock-capturing 
operator that locally reduces the order of accuracy near discontinuities. Three-dimensional viscous flow 
over the canopy of a Hermes computed by C'lialot et al. [32] is presented in Figures 3(a-b). The flow 
conditions corresponding to reentry are altitude of 60 km, = 20, angle of attack cv = 30° and Re/m 
of 120,890. The laminar solution was computed using equilibrium real gas hypothesis by the SUPG 
finite element method by an iterative implicit method that employs a block- diagonal preconditioned 
GMRES procedure. The mesh contains nearly 1 million tetrahedra. Figure 3a displays the surface grid 
and Figure 3b depicts the skin friction lines, illustrating the complex flow structure. 


5 Turbulence modeling 

In order to compute practical flows at large Reynolds numbers, turbulence has to be modeled. Some of 
the most popular turbulence models in use in structured grids such as the Baldwin-Lomax [10] become 
quite difficult to implement on unstructured grids because of their nonlocal nature, although this has 
been done [147, 104]. The trend is away from algebraic models and towards simple one and two equation 
field models, k - t turbulence models have been used with unstructured grids either by integrating the 
model to the wall [168] or by using a two-layer model [116, 32], where the two equation k - e model is 
replaced by a one equation model near the wall region. Marcum and Agarwal [100] have implemented 
two versions of k — c models in an unstructured grid code to compute turbulent flows over an ONERA 
M6 wing. They tested a low Reynolds number model that is integrated to the wall and a high Reynolds 
number model that is combined with wall- functions. Two fairly new one-equation models have become 
very popular, particularly for unstructured grid applications. These are the Baldwin-Barth [9] and the 
Spalart-Allmaras [159] models. The Baldwin-Barth model is derived from a simplified form of the k — e 
turbulence model. The Spalart-Allmaras model, on the other hand, is derived “from scratch” through 
dimensional analysis and intuition. Both the models have been tested extensively and are being used 
routinely. The cell-vertex finite volume scheme with a Galerkin procedure for computing the viscous 
terms has been used to compute viscous flows about high-lift configurations in conjunction with these 
turbulence models [6, 168]. In Figure 4, the surface pressure distributions and velocity profiles at three 
stations on the flap are displayed. The results are taken from the paper of Anderson and Bonliaus [6]. 
They use an implicit cell- vertex MUSCL scheme with the Spalart-Allmaras turbulence model, where the 
linear system at each time step is solved by a Gauss Seidel relaxation method. The results are compared 
with experimental data at two Reynolds numbers 5 x 10 6 and 9 x 10 6 with = 0.2 and angle of 
attack, a = 16°. The surface pressure profiles show' excellent agreement with the experimental data, 
while the velocity profiles show' good agreement except for the w'ake regions far downstream, where the 
grid resolution is inadequate. These and other quantitative comparisons with experimental data [168] 
indicate the level of maturity of unstructured grid technology. 


6 Time Discretization and Solution Strategies 


6.1 Steady-state solution techniques 


After discretizing Eqn. 
obtained: 


(1) in space, the following system of coupled ordinary differential equations is 


d(VMW) 


+ R(W ) = 0. 


dt 


( 6 ) 



Here IT is the vector of unknowns over the mesh points for a vertex-based formulation and over the 
cells for a cell-based formulation. V is the volume of the polyhedral control volume associated with 
the vertex/cell. In the case of a cell- vertex scheme storing pointwise quantities at the vertices, M is 
the mass matrix which represents the relationship between the average value in a control volume and 
the values at the vertices (the vertex representing the control volume and its nearest neighbors). It is 
only a function of the mesh and hence, a constant matrix for a static mesh. If a steady state solution 
is sought, time-accuracy is not an issue and M can be replaced by the identity matrix. This technique 
known as mass-lumping yields the following system of ordinary differential equations for the vector of 

unknowns IT : , TT - 

v^L + r(w) = 0 . (") 

at 

Cell-centered schemes (up to second order accuracy) and schemes dealing strictly with cell-averages (to 
any order of accuracy) do not yield a mass matrix and thus lead to Eqn. (7) for steady and unsteady 
problems. 

Eqn. 7 may be solved explicitly with linear multi-step methods of the forms described in [67, 
170]. The coefficients for these Runge-Kutta schemes are derived by considering a model problem and 
optimizing the coefficients to yield a large CFL number and good damping properties. Local time 
stepping, enthalpv damping (for Euler equations) and residual averaging are employed to accelerate 
convergence [74]. Even with this methodology, the convergence to steady state is usually unacceptably 
slow. Therefore, either multigrid methods or implicit schemes are required to accelerate the convergence. 

6.1.1 Multigrid methods 

The multigrid method has been demonstrated as an efficient means for obtaining steady -state solutions 
to the compressible Euler and Navier-Stokes equations on unstructured meshes in two and three dimen- 
sions. In this approach, convergence acceleration is achieved by time-stepping on successively coarser 
meshes. The principle behind this algorithm is that the errors associated with the high frequencies 
are damped by a carefully chosen smoother (e.g. a multi-stage Runge-Kutta scheme) while the errors 
associated with the low frequencies are damped on the coarser grids where these frequencies manifest 
themselves as high frequencies. In the case of structured grids, coarse grids are easily derived from a 
given fine grid by omitting alternate grid lines in each coordinate direction. In the case of unstructured 
grids, three different approaches can be adopted. 

The first approach begins with a coarse mesh definition and generates finer grids by refinement 
[131, 35]. The advantage is that the inter-grid operators become simple because of the nesting of grids. 
Another advantage is that this set-up can be utilized to advantage in an adaptive procedure, where the 
fine meshes are formed by adaptively refining the coarse meshes [35, 127]. The principal disadvantage of 
this approach is the dependence of the fine grid distribution on the coarse levels. The second approach 
uses non-nested unstructured grids either with a subset of fine grid points comprising the coarse grids 
or with completely independent coarse and fine meshes. This has been shown to be successful in both 
for inviscid and viscous flow computations [105, 101, 129, 24]. Both the approaches outlined above 
share a common problem, that of generating coarse grid levels. For complex geometries, especially in 
three dimensions, generating coarse grids that faithfully represent the complex geometries can become 
a difficult proposition. The requirement to generate not just one grid but multiple grids that preserve 
the geometry places too much of a burden on a user. The third approach circumvents this problem 
by generating coarse grids through agglomeration or fusing of fine grid control volumes, resulting in 
polyhedral coarse grid control volumes. This method was developed in [87] for cell-vertex schemes and 
in [158] for cell-centered schemes. This method has been further refined to deal with inviscid and viscous 
flows past complex configurations in both two and three dimensions [85, 177, 108, 109]. In Figures 5 (a-e) 
the results from using the agglomeration multigrid strategy are presented. Figure 5a depicts the surface 
grid employed to compute inviscid flow about a low-wing transport configuration [17/]. The mesh 
contains 804,056 vortices and approximately 4.5 million tetrahedra.. Figure 5b displays the computed 
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surface Mach contours for transonic inviscid flow using a seven-level agglomeration multigrid strategy. 
The freestream conditions for this case are = 0.77 and 1.116° incidence. Figure 5c shows the 
surface mesh employed for computing turbulent viscous flow over a partial span-flap wing experimental 
configuration [109]. The fine grid contains nearly 2.3 million vertices and nearly 13.6 million tetrahedra. 
The freestream conditions for this case are M = 0.2, incidence of 10°, and Reynolds number of 2 
million, corresponding to an approach condition. The solution obtained with a six-level agglomeration 
multigrid strategy is shown qualitatively in Figure 5d in the form of Mach contours on the wind-tunnel 
wall and density contours on the wing surface. Figure 5e shows the convergence histories for these two 
cases. The deterioration in convergence of the Navier-Stokes case is mainly due to the use of stretched 
meshes and indicates that the efficiency of Navier-Stokes solvers is far from satisfactory. 


6.1.2 Implicit schemes 

Implicit schemes for the compressible Euler and Navier-Stokes equations have been developed in order 
to accelerate the convergence to steady state. If the time derivative in Eqn. (7) is replaced by: 


dW W n+1 - W n 
dt ~ A t 


an explicit scheme is obtained by evaluating R(W) at time level n. An implicit scheme is obtained by 
evaluating R(W) at level n + 1. In the latter case, linearizing R about time level n, we obtain 


V OR \ 

Tt + m ' ■ ~ " 


AWi = (iy n+1 - w n y t . 


0 ) 

( 10 ) 


Eqn. (9) represents a large nonsymmetric linear system of equations for the updates of the vector 
of unknowns and needs to be solved at each time step. As At tends to infinity, the method reduces 
to Newton’s method. Direct solvers have been used solve Eqn. (9) yielding quadratic convergence, 
but these entail prohibitive costs and memory requirements and are impractical for three dimensional 
applications [176, 157]. Typically, due to storage considerations, only a lower order representation of the 
left-hand side is employed and the system is solved by iterative means. As a result of this approximation 
Eqn. (9) can never approach Newton’s method (with its associated quadratic convergence property) 
due to the mismatch of the right and left hand side operators. 

Thareja et al. [163] and Ha-ss an et al. [58] have utilized point-implicit iterative procedures. Fezoui 
[45], Batina [17], Anderson et al. [5, 7] and Slack et al. [157] have used a Gauss-Seidel relaxation 
technique. It is also possible to use more sophisticated techniques for the solution of the linear system, 
such as GMRES [149] and QMR [46]. Preconditioning of the matrices is critical to achieving good 
convergence with these methods. Some of the typical preconditioners are block-diagonal, symmetric 
successive over-relaxation (SSOR) and Incomplete LU factorization (ILU) preconditioners [111]. GM- 
RES with diagonal preconditioning has been used by Shakib et al. [151] to solve the linear systems 
arising out of a finite element discretization of the Euler equations. Slack et al. [157] and Whitaker [184] 
have also used GMRES with diagonal preconditioning in two and three-dimensional applications. Slack 
et al. [157] have observed when solving the two-dimensional Euler equations that the preconditioned it- 
erative methods perform better than the other methods as the number of elements in the mesh increases. 
Venkatakrishnan and Mavriplis [178] tested a family of implicit schemes for solving the two-dimensional 
compressible Navier-Stokes equations on unstructured meshes. They concluded that GMRES with ILU 
preconditioning (GMRES/ILU) was superior to other implicit schemes over a range of problems sizes 
and flow conditions and was competitive in terms of CPU time with multigrid methods. Luo et al. 
[99] used GMRES/ILU to compute two- and three-dimensional flows. A drawback of all the methods 
based on linear algebra is that the memory requirements become severe which seriously limits the sizes 
of the problems that can be solved, especially in three dimensions. It is possible to implement GMRES 
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in a matrix- free form [29] where the matrix- vector product is replaced by a' finite difference expression 
involving residual evaluations. Another advantage of the matrix-free approach is that the higher order 
discretization can be employed for the left-hand side and the matrix- vector product [30]. The drawback 
of this approach is that powerful preconditioners like SSOR and ILT cannot be used since they require 
the matrix to be explicitly available. A compromise is to settle for a matrix-free form of GMRES with 
diagonal preconditioning as done in [75]. 

Another implicit method that has been investigated is based on the use of linelets [59]. The key idea 
is to cover the domain with a set of lines and the scheme is made implicit along these lines. Thus tri- 
diagonal systems need to be solved which can be accomplished efficiently. A drawback of this method 
is that the convergence is sensitive to the orientation of the linelets with respect to the direction of 
strongest coupling. 

6.2 Unsteady Problems 

While the solution techniques for computing steady flows have evolved to a high degree, those for dealing 
with unsteady flows have lagged behind. This is partly due to fact that prediction of aerodynamic 
properties in steady transonic flow (cruise conditions) has been the driving force, with the unsteady 
effects being important only in extreme conditions, such as flutter and stall. It is anticipated, however, 
that unsteady time-accurate simulations will provide the next great challenge for CFD. 

Recall that after discretization in space, a system of coupled ODE’s results (Eqn (6)). If we assume 
that the system is decoupled i.e. M = /, the simplest way to solve the system of ODE s is to use an 
explicit scheme. Low-storage Runge-Kutta methods [67, 170] and various predictor-corrector schemes 
fall under the category of explicit schemes. These schemes typically require only simple updates. They 
place restrictions on time step due to the CFL condition, which become severe in the case of semi- 
discrete schemes satisfying monotonicity principles. However, explicit schemes may be the schemes of 
choice for certain unsteady applications when the time scales of interest are small or more precisely, 
that they are comparable to the spatial scales. The grid should be clustered only in regions of interest; 
otherwise, the size of the explicit time step could be unnecessarily small. The situation can be improved 
by the use of time-step sequencing [83, 138], where different cells take varying number of local time 
steps to get to a particular time level. 

When the mass matrix is present, even the explicit schemes require a matrix inversion. Often, either 
the mass-matrix is lumped for convenience or at best a few Jacobi iterations are carried out [124, 36]. 
Neither approach is entirely satisfactory, especially if large time steps, permitted by multi-stage explicit 
schemes and implicit schemes, are used. Donea [41] realized that the presence of the mass matrix in 
his formulations will make it unattractive and proposed a simple two-pass procedure. While he proved 
that the formal order of accuracy remained the same, he also showed that the dissipation and dispersion 
properties were somewhat compromised. Recently, it has been observed [179] that the mass matrix can 
be lumped without adverse consequences if the cell- vertex scheme possesses only second or lower order 
accuracy. 

When an implicit scheme is used to solve for unsteady flows, the linear system Eqn. (9) is modified 
by replacing the diagonal matrix V by the matrix MV where M is the mass matrix and V is the volume 
of the ceil. However, it is not enough to solve this linear system. One has to drive the unsteady residual 
to zero or at least to truncation error. This is usually done by employing inner iterations [139, 136]. It 
is the role of these inner iterations to eliminate errors due to the factorization (if any is carried out), 
linearization, and errors arising from employing a lower order approximation on the implicit side. The 
number of inner iterations required may be large depending on the flow situation and the size of the 

time step employed. 

Brennis and Eberle [27] and Jameson [68] have advocated a different approach for deriving an 
efficient implicit scheme for unsteady flows. The idea is to define an unsteady residual, following a 
backward difference approximation to the time derivative and then use either a relaxation strategy or a 
multigrid technique to drive the unsteady residual to zero. The significant advantage of this approach 
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when multigrid is used to solve the nonlinear problem is that it incurs no storage overheads associated 
with traditional implicit schemes, and is particularly attractive for unstructured grid computations in 
three dimensions. It allows the time step to be determined solely based on flow physics. This method 
has been used to compute two- and three-dimensional flows over airfoils and wings [68, 112, 4] using 
structured grids. Vassberg [172] has applied this method to compute flow solutions over oscillating 
airfoils using unstructured grids. This approach has also been used in [179] in conjunction with the 
agglomeration multigrid strategy to solve two-dimensional unsteady flows where the inversion of the 
mass matrix is accomplished indirectly during the multigrid procedure. 

7 Adaptive grids 

Solution adapted grids are increasingly being used to compute complex unsteady and steady flows. 
There are three distinct ways the grid can be adapted to the solution. These are r-refinement, h- 
refinement and p-refinement. In r-refinement, the nodes are simply redistributed so that regions of 
importance are better resolved. In h- refinement or mesh-enrichment, the cells are locally subdivided 
or merged or in some instances, a complete remeshing is done to reduce the grid spacing in regions of 
interest. In p-refinement, the degree of the basis function is adjusted locally to match the variation 
in solution. For a survey on adaptive mesh refinement techniques, the reader is referred to the review 
article by Powell et al. [135]. 

R-refinement is probably the simplest in concept, but is burdened with practical difficulties in multi- 
dimensions especially when dealing with highly stretched grids. For a survey of mesh point movement 
methods, see the review article by Haw ken et al. [60]. The difficulties include skewness, crossing of 
lines, arbitrarily small cell volumes etc. Some progress has been made in dealing with these issues 
at least for inviscid flows, where the cell-aspect ratios are not too severe [122]. The advantage of r- 
refinement is that if a valid grid results from it, all that is required is the interpolation of variables from 
the old to the new grid. This could be done in a conservative manner if desired. A way to avoid the 
interpolation, which introduces errors that could accumulate, is to introduce the grid movement terms 
in the governing equations Eqn. (1). These terms need to be discretized carefully so that freestream 
is preserved. The Geometric Conservation Law (GCL) [165, 186] formalizes this procedure. Recently, 
r-refinement has been used to great advantage with Roe's upwind scheme [143] to obtain “fitted” shock 
resolution by Paraschivoiu et al. [123], Parpia and Parikh [126], and van Rosendale [171] by aligning 
edges of the triangulation with discontinuities. The “fitting” is done in a shock-capturing framework 
by utilizing that property of Roe’s scheme which allows isolated discontinuities aligned with the mesh 
to be captured exactly. 

Il-refinement is by far the most popular means of adaptation in compressible flow T s. This is especially 
true for inviscid flows dominated by interactions of shock waves where p-refinement techniques are of 
limited value. In the Adaptive Mesh Refinement (AMR) approach, h-refinement is employed with body- 
fitted structured grids to adapt to the flow solution [20, 138]. Instead of body-fitted grids, Cartesian 
meshes can be used to adapt to complex geometries in addition to adapting to flow features [21, 37, 137]. 
Since some of the approaches, e.g. [37, 152] do not make use of any structure and employ indirect 
addressing, these methods may be classified as unstructured. The AMR framework has been used 
to compute many complex inviscid unsteady flows that, exhibit disparate spatial and temporal scales. 
Unstructured grids, composed of triangles or tetrahedra, are also particularly suited for h-refinement. In 
the case of both AMR and unstructured grids, the regions of interest are first identified either through 
a combination of heuristic criteria such as density gradients (undivided) or through estimation of the 
truncation error. In the case of AMR, typically cells are subdivided in some fixed ratio, and these are 
sometimes grouped into logically rectangular patches. In the case of unstructured grids, the regridding 
procedure is more involved. If a Delaunay triangulation is used as the initial grid, one can use Watson’s 
incremental algorithm [183] to obtain a new Delaunay triangulation. Other procedures such as edge- 
swapping may also be used to restore the Delaunay triangulation. In the case of other triangulations, 


11 



cells are subdivided and smoothed using a few passes of a Laplacian-type smoother. Peraire et al. [130] 
advocate the regeneration of the grid by the advancing front method for steady state applications after 
identifying regions of interest by error estimation. In the case of nodal schemes in three dimensions, 
when edges are tagged and the cells are to be refined or derefined based on these tagged edges, the 
possibilities are numerous [79]. It helps to use some rules to limit the cases [35]. Several investigators 
have employed multigrid procedure for steady state computations using the grids generated by solution 
adaptation [102, 127, 35, 26]. Figure 6a shows the initial surface grid for an isolated engine inlet. 
The initial grid contained 3536 nodes. After three adaptive refinement cycles, the final mesh shown in 
Figure 6b was obtained. The final mesh contained 61,002 nodes. This computation done by Connell 
and Holmes [35] utilized the adaptive multigrid procedure. When adapting grids to flow- features, a 
problem that, may occur is that certain regions are overly resolved at the expense of other regions, 
resulting in seemingly sharp but incorrect solutions. Warren et al. [182] have proposed and tested a 
simple modification to the error indicators that alleviates this problem. 

In the case of transient problems, adaptation is performed much more frequently and therefore the 
regridding process is required to be efficient. While this is easily accomplished efficiently by subdivision 
of cells in the AMR framework, the problem is more involved in the case of unstructured grids [94, 141]. 
The problem of flow past bodies in relative motion has also been addressed in the literature [93, 167, 
57, 82, 31]. Tvpically, mesh point movement and efficient mesh restructuring are employed to obtain 
valid, good-quality grids about the moving bodies. Figures 7(a-c) show the results from the simulation 
of store separation from F-117. The freestream conditions are = 0.84 and incidence of 2.73°. 
The projectiles are forcibly ejected at different velocities. This multi-body problem was computed by- 
Baum et al. using the Arbitrary Euler Lagrangian adaptive unstructured methodology developed in 
[18]- Because of symmetry, the motion of only two of the projectiles is simulated. Figures 7a and 7b 
show the Mach contours from two different perspectives and Figure 7c displays the pressure contours 
at a particular instant of motion. The complex shock patterns emanating from the projectiles and their 
signatures on the plane of sy r mmetry r may be observed in Figure 7c. 


8 Higher order methods 

The use of methods possessing higher than second order accuracy for solving the compressible Navier- 
Stokes equations on unstructured grids is not yet commonplace and is a topic of active research. The 
finite element method achieves higher order accuracy by the use of higher order elements. Bey and 
Oden [22] have used a discontinuous Galerkin method to obtain up to 4th order accuracy for smooth 
flows using structured adaptive grids. In the finite volume community, Barth and Fredrickson [14] 
derived general conditions for a scheme to be higher order accurate that involve the reconstruction 
of variables satisfying the properties of conservation of mean, ^-exactness and compact support, k- 
exactness refers to the property’ of being able to reconstruct exactly polymomials of degree ^ A. The 
key idea is to extend the support of the scheme to enable the coefficients in the polymomial to be 
determined. They also proposed a minimum-energy reconstruction that used a larger support and 
utilized a least-squares method to evaluate the coefficients. By combining this reconstruction procedure 
with higher order quadrature, they demonstrated the effectiveness of both h- and p-refinements in 
reducing the errors in smooth flows. More recently, Chakravarthy et al. [31] have proposed a similar 
approach to achieving higher order accuracy-. Equivalent to requiring A-exactness, they propose that 
the mean be conserved in a neighborhood of the point in question. Mitchell [115] has proposed a quasi- 
quadratic reconstruction procedure that appears to ydeld better results than those obtained w-ith linear 
reconstruction. Harten and Chakravarthy [55] have proposed a framework for applying Essentially 
Nonoscillatory (ENO) schemes [56] on unstructured meshes. They present two techniques for adapting 
the stencils during the reconstruction procedure, which also give preference to a centered stencil in 
smooth regions. The first technique considers several candidate- stencils and picks the one that represents 
the function the smoothest. The second technique is a hierarhical one that incrementally augments the 
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stencil as derivatives are required during the reconstruction procedure. Abgrall [2] has also proposed a 
hierarchical ENO scheme and has tested it by computing a variety of inviscid compressible flows. 

A problem with the reconstruction methods described above is that since the support grows with the 
order of the scheme, computational cost and boundary treatment become significant. Halt and Agarwal 
[53] have proposed and tested two methods for realizing higher order accuracy while keeping the support 
compact. These techniques are based on deriving equations for either the derivatives or the moments 
of the governing equations. They demonstrate up to 4th order accuracy with both these methods for 
smooth flows. They conclude that the method based on moments is more suitable for dealing with shocks 
and boundary conditions. Barth [13] has also proposed a different quadratic reconstruction procedure. 
Similar to a finite element procedure, additional degrees of freedom are placed at the mid-points of edges 
in a triangulation. The gradient and Hessian required for quadratic reconstruction are again obtained 
by a least-squares procedure. Figures 8(a-d), provided by Barth [13], illustrate the results obtained with 
various discretizations for the circular advection problem: 


u t + (yu) x - ( xu) y = 0. (11) 

Discontinuous initial data is specified along an interior cut line. The exact solution is a solid body 
rotation of the cut-line data. The triangular grid is displayed in Figure 8a. Solutions obtained with 
the finite volume scheme employing piecewise-constant, piecewise- linear and piecewise-quadratic recon- 
structions are depicted in Figures 8 (b-d). The significant improvement with higher order reconstruction 
is apparent. Barth [13] also states that quadratic reconstruction is about 7 times as expensive as linear 
reconstruction; of this a factor of 4 arises from the introduction of grid points at the mid-points of the 
edges. This “h- refinement” contributes in part to the improvement realized with quadratic reconstruc- 
tion. 

9 Hybrid discretizations 

While triangulations in two and three dimensions offer flexibility in adapting to complex geometries, 
it is not clear if they are needed near the wall regions. Steger [161] when introducing the “thin- 
layer” approximation in 1978 stated, “One does not generally have sufficient computer to resolve the 
viscous terms except in a thin-layer near the body”. This statement holds true even today for practical 
aerodynamic computations. Thus highly stretched triangulations are required in the near-wall regions. 
They can be generated by employing a mapping procedure [103] where a Delaunay triangulation is 
performed in the mapped plane which is then mapped back to the physical space to yield the desired 
stretched triangulation. This method has also been used for adaptation in viscous flows [181]. Aftosmis 
et al. [3] have addressed the important issue of element shapes by examining the accuracy obtained 
on test problems by using various triangular element shapes and quadrilaterals. They conclude that 
using stretched triangulations near the body does not yield any improvement in solution accuracy 
over using quadrilaterals, and is therefore inefficient because of the extra edges in the triangulation. 
They recommend a method that removes unnecessary edges (diagonals) from boundary layer regions. 
Rather than adopt this two-step procedure, where the triangular grid is first generated followed by the 
elimination of unnecessary edges, the hybrid method eliminates the need for stretched triangulations 
all together. This method makes use of a structured or a semi-structured grid in the near- wall regions. 
In two dimensions, a structured body-fitted grid is used near the body whereas in three dimensions 
a prismatic grid is used since the surface grid is triangular. Nakahashi et al. [117, 118], Kaflinderis 
et al. [80, 128], Connell and Braaten [34] have advocated this approach. The generation of prismatic 
elements requires careful marching away from the surface so that lines do not cross. Melton et al. [113] 
have used a prismatic grid with a background Cartesian mesh to compute three-dimensional flows. The 
structure afforded could be effectively used to tailor line-implicit algorithms, although this requires that 
the various elements be processed separately. Ideally, one data structure (edge-based) should be able to 
process all elements regardless of their shapes for an automatic procedure. Figure 9a shows the surface 
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grid for an adapted hybrid prismatic-tetrahedral grid about a hemisphere. The initial grid contained 
3000 triangular faces on the wall and 20 layers of prisms and 70,000 tetrahedra. The freestram conditions 
for supersonic laminar flow are M^o — 1-4 and Rc — 1000. After the dual adaptation procedure, where 
the prisms are directionally refined with the normal resolution unaltered and the tetrahedra are divided 
into either 2,4, or 8 tetrahedra, the final grid contains 7000 triangular faces and 20 layers of prisms and 
225,000 tetrahedra. The signature of the directional refinement of the prisms may be observed on the 
hemispherical surface. Figure 9b displays the surface Mach contours on the symmetry plane for this 
flow. These figures have been provided by Parthasarathy et al. [128]. The explicit solver also makes 
use of an adaptive multigrid procedure. 


10 Parallel computing 

Computational fluid dynamics as its name implies is inevitably linked to computing issues. Among 
these are processing power, memory technology, networking and accessibility'. Ability to compute the 
solutions to problems in finite time always being the goal, CFD has benefited immensely from the rev- 
olution that has taken place in the last 15 years in these areas. Vector supercomputers have provided 
much of the computing power that has been harnessed to compute complex three-dimensional flow's. 
It is anticipated that distributed-memory parallel computers will offer the next cost-effective leap for- 
ward in terms of computing pow r er. Much research has been carried out in the area of unstructured 
grid computations on parallel computers. Some of the issues that arise are partitioning of the grid, 
message patterns, data structures and design of parallel algorithms. Partitioning of unstructured grids 
for parallel computing has been investigated by a number of researchers. The methods can be broadly 
classified into geometry-based [23, 114, 44] and graph-based [133, 156] algorithms. It appears that the 
graph-based algorithms, in particular, the spectral bisection technique of Pothen et al. [133], wdiile 
being computationally more expensive, yield much better partitions. Unstructured grid flow' solvers 
have been implemented on various machines [180, 107, 75, 43, 140]. These studies have shown that 
good performance may be obtained by paying careful attention to the issues listed above. Regarding 
the flow solver, explicit schemes and point Jacobi implicit schemes possess almost complete parallelism, 
except for communication at the inter-processor boundaries. Multigrid methods appear to have ade- 
quate parallelism for coarse-grained parallel computers; impressive performances have been reported by 
Mavriplis et al. [107]. Johan et al [75], Ramamurthi et al. [140] and Venkatakrishnan [175] have showm 
that implicit schemes can be carefully designed to yield good performance when solving unstructured 
grid problems on parallel computers. While a matrix-free GMR.ES w'ith nodal block preconditioner 
that has almost complete parallelism except for the communication at the inter-processor boundaries 
is used in [75], methods that are implicit within processors but explicit across processors are used in 
[140, 175]. Thus, as the number of processors increases, the latter methods exhibit a degradation in 
convergence. In [175] it is demonstrated that convergence rate of such a method can be improved by- 
using in addition a global coarse grid approximation, where each processor is represented by one grid 
point. In general, it appears that the best parallel algorithm is also the best serial one, with some 
possible loss in convergence in the case of implicit schemes in a multi-processor environment. 


11 Conclusions 

In this paper, unstructured grid flow' solvers for the compressible Navier-Stokes equations have been 
surveyed. Significant progress has been made in the areas of spatial and temporal discretization, adaptive 
and parallel algorithms. Newly developed field equation turbulence models seem to mesh nicely with 
the unstructured grid framework. Based on these developments, it is safe to say that unstructured 
grid technology is almost on par with structured grid technology, although encumbered with additional 
memory and computational costs. This overhead has to be balanced with the ability to compute flow's 
over complex geometries and the ease of adaptation. 
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The areas that require further research are thus the same as in structured grids. These include better 
and faster implicit/multigrid flow solvers that are insensitive to cell aspect ratios and grid stretching, 
improved higher order discretization techniques, better a priori or a posteriori error estimates and 
parallel algorithms. The status of unsteady flow solvers is far from satisfactory either for structured or 
unstructured grids. The use of hybrid meshes needs to be explored further to utilize data structures 
that can handle any type of tessellation. As unsteady flows and aeroacoustics are emerging as areas 
of interest, nondissipative schemes that also minimize dispersion have to be designed to be applicable 
to unstructured grids, which are nonuniform in general. In this respect as well as in the areas of error 
estimation and higher order schemes, there should be more interaction between finite volume and finite 
element communities. 
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Figure 1. Mach contours for inviscid transonic flow over Boeing 747-200 using an unstructured tetra- 
hedral grid (Mqo = 0.8, a = 2.73°) [72]. 
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Figure 3. Hypersonic flow over the canopy of a Hermes [32]. (a) Surface grid. (b). Skin friction lines 
on the surface. 
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■Fig. 5b. 
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Fig. 5e. 

Figure 5. Computations done using the agglomeration multigrid strategy [177, 109]. (a) Surface grid for 
the low wing transport (LWT) case, (b) Mach contours on the surface for inviscid transonic flow over 
the LWT. (c) Surface grid for the partial span-flap experimental configuration, (d) Mach contours on 
the wind-tunnel wall and density contours on the surface for turbulent flow, (e) Convergence histories 
for the inviscid and the viscous cases. 
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Figure 6. (a) Initial surface grid for an engine inlet, (b) Grid obtained through adaptive refinement 

[35]. 




Fig. 7c. 


35 

















REPORT DOCUMENTATION PAGE 


Form Approved 
OMB No 0704-0188 


PublicTe porting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, 
gathering and maintaining the data needed, and completing and reviewing the collection of information Send comments regarding this burden estimate or any other aspect of ilm 
collection of information, including suggestions for reducingthis burden, to Washington Headquarters Services. Directorate for l nforma ;'^pP e ^ i0 ^ 3 " d 9 J P ° n A‘rl\l n { tiltTSOn 
Davis Highway, Suite 1204, Arlington^VA 22202-4302, and to the Office of Management and Budget, Paperwork Redu ction Project (0704-0188), Washington, DC 20503 


1. AGENCY USE ONLY (Leave blank) 2. REPORT DATE 

February 1995 


4. TITLE AND SUBTITLE 

A PERSPECTIVE ON UNSTRUCTURED GRID FLOW SOLVERS 


3. REPORT TYPE AND DATES COVERED 

Contractor Report 


5. FUNDING NUMBERS 


C NA31-19480 
WU 505-90-52-01 


6. AUTHOR(S) 

V. Venkatakrishnan 


7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Institute lor Computer Applications in Science 
and Engineering 

Mail Stop 132C, NASA Langley Research Center 
Hampton, VA 23681-0001 


9. SPONSORING/MONITORING AGENCY NAME{S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, VA 23681-0001 


8. PERFORMING ORGANIZATION 
REPORT NUMBER 

ICASE Report No. 95-3 


10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

NASA CR-195025 
ICASE Report No. 95-3 


11. SUPPLEMENTARY NOTES 

Langley Technical Monitor: Dennis M. Bushnell 
Final Report 

Submitted to AIAA Journal 


12a. DISTRIBUTION/AVAILABILITY STATEMENT 

U nclassified-L^ nlimited 


12b. DISTRIBUTION CODE 


Subject Category 64 


13. ABSTRACT (Maximum 200 words) . , 

This survev paper assesses the st atus of compressible Euler and N avier-Stokes solvers on unstruct ured grids. JJiiterent 
spatial and temporal discretization options for steady and unsteady flows are discussed. The integration of these 
components into an overall framework to solve pract ical problems is addressed. Issues such as grid adaptation, higher 
order methods, hybrid discretizations and parallel computing are briefly discussed. Finally, some outstanding issues 
and future research directions are presented. 


14. SUBJECT TERMS 

Unstructured; Finite Volume; Finite Element; Adaptive Grids 


17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

NSN 7540-01-280-5500 


15. NUMBER OF PAGES 

39 


16. PRICE CODE 

A03 


118. SECURITY CLASSIFICATION! 19. SECURITY CLASSIFICATION] 20. LIMITATION 


OF THIS PAGE 

Unclassified 


OF ABSTRACT 


OF ABSTRACT 


Standard Form 298(Rev. 2-B9) 

Prescribed by ANSI Std 239- 1 8 
298-102 










National Aeronautics and 
Space Administration 
Langley Research Center 
Mail Code 1 80 
Hampton, VA 23681-00001 

Official Business 
Penalty for Private Use, S300 


BULK RATE 

POSTAGE & FEES PAID 

NASA 




